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Abstract 



Ph. 

For many complex dynamical systems, a separation of scales prevails between the 
Q i (microscopic) level of description of the available model, and the (macroscopic) level at 

which one would like to observe and analyze the system. For this type of problems, an 
"equation-free" framework has recently been proposed. Using appropriately initialized 
^ ' microscopic simulations, one can build a coarse-grained time-stepper to approximate a 

^^' time-stepper for the unavailable macroscopic model. Here, we show how one can use 

this coarse-grained time-stepper to compute coarse-grained traveling wave solutions of 
a lattice Boltzmann model. In a moving frame, emulated by performing a shift-back op- 
eration after the coarse-grained time-step, the traveling wave appears as a steady state, 
^ ' which is computed using an iterative method, such as Newton-GMRES. To acceler- 

.,—1. ' ate convergence of the GMRES procedure, a macroscopic model-based preconditioner 

is used, which is derived from the lattice Boltzmann model using a Chapman-Enskog 
expansion. We illustrate the approach on a lattice Boltzmann model for the Fisher 



^^3 ' equation and on a model for ionization waves. 
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1 Introduction 

There is an established algorithmic infrastructure to study the long-term dynamical fea- 
tures of systems of partial differential equations (PDEs), such as steady states or periodic 
solutions. When only a simulation code (a time-stepper) is available, algorithms such as 
the recursive projection method (RPM) [34] and Newton-Picard [24, 26] can locate steady 
states, as well as their stability, and perform a continuation for changing values of the param- 
eters. These methods project the Jacobian onto the eigenspace corresponding to the slowly 
decaying modes, which is typically low-dimensional. In this subspace, a Newton iteration is 
performed; in the orthogonal complement, Picard iterations (time-stepping) converge fast 
enough to the steady state. Alternatively, so-called Jacobian-free Newton-Krylov methods 
[20] solve the linear system for each Newton correction by means of an iterative method, 
such as GMRES, for which an appropriate preconditioner is crucial. Both standard pre- 
conditioning techniques, such as incomplete LU factorization (ILU) and multigrid, as well 
as application-specific physics-based preconditioners have been proposed, see [20] for an 
overview and references. If required, the stability can be computed as a post-processing 
step [21]. As a common feature, all these methods only use selected matrix- vector products 
with the system's Jacobian, which are estimated using the time-stepper with several nearby 
initial conditions. 

Unfortunately, a low-dimensional macroscopic PDE is often not able to capture all 
detailed physical interactions accurately. In such cases, one needs to resort to a more 
microscopic description. For instance, the dynamics of a system of colliding particles with 
interactions that depend sensitively on the relative particle velocities can, in general, not 
be modeled by a reaction-diffusion equation for the particle density. One example, which 
forms the main motivation for the present paper, is the impact ionization reaction, where 
each collision of an electron with a neutral atom or molecule creates an additional electron 
when the relative velocity is above a certain treshold. Such a dynamical system should be 
modeled through a phase space evolution law, e.g. a Boltzmann equation. 

Nevertheless, a clear separation in time-scales is often present in the microscopic model: 
on fast time-scales, the microscopic variables equilibrate with respect to a few macroscopic 
variables, while these macroscopic variables themselves evolve on much slower time-scales. 
When this is the case, a macroscopic model should conceptually exist. However, it might 
be difficult (or impossible) to derive a closed expression from the underlying microscopic 
model without introducing assumptions that are hard to justify. 

For such models, there is an active current interest in the development of so-called 
equation-free methods to study the long-term behavior [19]. The key idea, which was first 
proposed in [38], is to construct a coarse-grained time-stepper for the unavailable macro- 
scopic equation as a three step procedure: (1) lifting, i.e. the creation of appropriate initial 
conditions for the microscopic model, conditioned upon the macroscopic state at time t*; (2) 



simulation, using the microscopic model, over the time interval [t*,t* + 6t]; (3) restriction, 
i.e. the extraction of the macroscopic state at time t* + 6t. The result is a coarse time-dt 
map, which can be used to estimate the matrix- vector products that are required in an 
RPM, Newton-Picard or Newton-Krylov method. 

Based on RPM, coarse-grained bifurcation analysis has already been used in a number of 
applications [17, 36], and also allows to perform other system-level tasks, such as control and 
optimization [35]. In this paper, we will investigate the use of Jacobian-free Newton-Krylov 
techniques on a model problem concerning traveling wave solutions of lattice Boltzmann 
models. Traveling waves are solutions that move with constant speed without changing 
shape; in a co-moving frame, they appear as steady state solutions. We construct a coarse- 
grained time-stepper in this co-moving frame by performing a shift-back operation after each 
coarse-grained time-step, and compute its fixed points using a Newton-GMRES procedure. 
To accelerate convergence of the GMRES procedure, we build a preconditioner based on 
an approximate PDE model, which is derived from the lattice Boltzmann model through 
a Chapman-Enskog expansion. As a consequence, the method described here could more 
appropriately be called equation- assisted, rather than equation-free. We expect that the 
techniques described here can be applied in other applications where particle based methods 
are necessary to describe the dynamics. 

This paper is organized as follows. In section 2, we briefly review the basic properties 
of the coarse-grained time-stepper. Subsequently, we outline the model problems that will 
be used throughout the text in section 3. The non- linear system, of which the traveling 
waves are the solution, is constructed in section 4, and the preconditioned Newton-GMRES 
method is discussed in section 5. Section 6 contains a detailed numerical study of the 
convergence properties of the method. Finally, we conclude in section 7, which contains a 
discussion of the computational complexity and some final remarks. 

2 Coarse-grained time-stepper 

We briefly review the coarse-grained time-stepper, as it was introduced by Kevrekidis et al. 
[19]. To this end, we consider an abstract microscopic evolution law, 

dtu{^,t) = f{u{^,t)), (1) 

in which w(x, t) represents the microscopic state variables, x G Dm and t are the microscopic 
independent variables, and dt denotes the time derivative. We assume that a macroscopic 
model, denoted by 

5it/(X,t)=F(C/(X,t)), (2) 

conceptually exists, but is unavailable in closed form. In equation (2), f/(X, t) represents the 
macroscopic state variables, and X G Dm and t are the macroscopic independent variables. 



We introduce a time-stepper s for the microscopic evolution law (1), 

■u(x,t + dt) = s(n(x,t);dt), (3) 

where dt is the size of the microscopic time-step, and the aim is to obtain a coarse-grained 
time-stepper S for the variables t/(X, t) as 

U{X,t + 6t) = S{U{'X,t);6t), (4) 

where 5t denotes the size of the coarse-grained time-step, and the bars have been intro- 
duced to emphasize the fact that the time-stepper for the macroscopic variables is only an 
approximation of a time-stepper for (2), since this equation is not explicitly known. 

To define a coarse-grained time-stepper (4), we need to introduce two operators that 
make the transition between microscopic and macroscopic variables. We define a lifting 
operator, 

fi : C/(X, t) ^ n(x, t) = ;u(C/(X, t)), (5) 

which maps macroscopic to microscopic variables, and its complement, the restriction op- 
erator, 

M:u{j^,t)^U{X,t)=M{u{^,t)). (6) 

The restriction operator can often be determined as soon as the macroscopic variables are 
known. For instance, when the microscopic model consists of an evolving ensemble of many 
particles, the restriction typically consists of the computation of the first few low order 
moments of the distribution (density, momentum, energy), which are considered as the 
appropriate macroscopic variables C/(X, t), in terms of which a closed macroscopic equation 
can be written. The assumption that a macroscopic equation closes at the level of these low 
order moments, implies that the higher order moments become functionals of the low order 
moments on time-scales which are fast compared to the overall system evolution (slaving). 

The construction of the lifting operator is usually more involved. Again taking the 
example of a particle model, we need to define a mapping from a few low order moments 
to initial conditions for each of the particles. We know that the higher order moments 
of the resulting particle distribution should be functionals of the low order moments, but 
unfortunately, these functionals are unknown (since the macroscopic evolution law is also 
unknown). Several approaches have been suggested to address this problem. One could for 
instance initialize the higher order moments randomly. This introduces a lifting error, and 
one then relies on the separation of time-scales to ensure that the higher order moments relax 
quickly to a functional of the low order ones (healing) [14, 27, 36] (see also [30, 39]). We note 
that, in some cases, this approach produces inaccurate results [40]. In fact, to initialize the 
higher order moments correctly, one should perform a simulation of the microscopic system 
subject to the constraint that the low order moments are kept fixed. How this can be done 
using only a time-stepper for the original microscopic system, is explained and analyzed in 
[12, 13, 41]. We will briefiy discuss the lifting step for our model problems in section 4.1. 



Given an initial condition for the macroscopic variables U(X.,t*) at some time t*, we 
can construct the time-stepper (4) in the following way: 

1. Lifting. Using the lifting operator (5), create appropriate initial conditions u{x,t*) 
for the microscopic time-stepper (3), consistent with the macroscopic variables. 

2. Simulation. Use the microscopic time-stepper (3) to compute the microscopic state 
■u(x,t) forte [t*,t* + 6t]. 

3. Restriction. Obtain the macroscopic state C/(X, t* + 6t) from the microscopic state 
u(x,t* + 6t) using the restriction operator (6). 

Assuming 5t = kdt, this can be written as 

C7(X, t + 6t) = S(i7(X, t),6t) = M{s''{fi{U{X, t)),dt)), (7) 

where we have represented the k microscopic time-steps by a superscript on s. If the 
microscopic model is stochastic, one may need to perform multiple replica simulations, 
using an ensemble of microscopic initial conditions, to obtain an accurate result with a 
sufficiently low variance. 

3 Model problems 

We now briefly discuss the origins of the lattice Boltzmann method and the relation with 
the Boltzmann equation; we also introduce our model problems. 

3.1 Microscopic and macroscopic model 

We consider systems that, on a molecular level, consist of particles whose position and ve- 
locity are governed by two processes: free flight and collisions. We introduce the probability 
f {x , V , t)dxdv , which represents the fraction of particles with position and velocity in the 
infinitesimal domain [x, x + dx] x [v,v + dv] at time t. The evolution of / is governed by a 
so-called kinetic equation, 

dtf{x, V, t) + vd^fix, V, t) + Fd,f{x, V, t) = QU, /), (8) 

where Q(/, /) is the collision integral, and F is an external force term. The first equation of 
this type was the Boltzmann equation for moderately rarefied gas flows with F = [2]. Of 
course, when multiple species are present, (8) becomes a system of equations. Throughout 
this paper, we confine ourselves to problems in one space dimension. 

Usually, one introduces a kinetic model for Q{f,f) to simplify equation (8) [15]. A 
standard choice is the non-linear Bhatnagar-Gross-Krook model (BGK) [1], 

dtfix, V, t) + vdj{x, V, t) + FdJix, V, t) = _ /(^iM:i/!!(^iM. (9) 



Here, the collisions are interpreted as a relaxation to local equilibrium distribution, with a 
characteristic relaxation time r. The choice of the equilibrium distribution f^'^{x, v, t) (and 
therefore of the collision integral Q{f,f)) determines the physics of the system. 

Although the kinetic equation (8) and its BGK-approximation are able to describe the 
evolution of a wide range of physical systems, a numerical simulation quickly becomes in- 
tractable because of the high dimensionality. However, one can often obtain an approximate 
macroscopic description. We define the moments of the distribution function as 

/oo 
v'f{x,v,t)dv, i = 0,1,2,..., (10) 

-oo 

where m is the particle mass. The lowest order moments correspond to the density (z = 0), 
momentum {i = 1) and energy (i = 2), which we write as 

p(x,t)=/>W(x,t), cl)ix,t) = p^^\x,t), i{x,t) = p^^\x,t)/2. (11) 

Using the assumption that the deviation from local equilibrium is sufficiently small, one 
performs an asymptotic expansion (the Chapman-Enskog expansion, [22]) to obtain a closed 
description for the evolution of a number of low order moments. If we define the macroscopic 
variables as U{x,t) = {/)'*' (x, t)} ._, we obtain an approximate PDE of the form 

dtU{x, t)=F (U{x, t),d^U{x, t), . . . , diU{x, t)) , (12) 

which depends on the first d spatial derivatives. Equation (12) is a good approximation of 
the system dynamics when there is a fast decay of the higher order terms in the Chapman- 
Enskog expansion. 

Consider as an example the advection-diffusion equation, 

dtp{x, t) + cdxp{x, t) = dx{Dda:p{x, t)), (13) 

with transport coefficient c and diffusion coefficient D. This equation can be derived from 
equation (9) with F = hy defining the equilibrium distribution as 

nx, V, t) = p{x, t)^exp{-Xiv - cf), 

with A = m/2kT, where T is temperature and k is the Boltzmann constant, and adding the 
conservation constraint 

/•oo 

{f{x,v,t) - g{x,v,t))dv = 0. 



A straightforward derivation reveals that the macroscopic behaviour of (9) is described by 
equation (13) when we choose r = 2DX [44]. Different equilibrium distributions can be 
used to obtain the Burgers' equation, the Euler equations, the Navier-Stokes equations, 
etc. One can also add chemical reactions in the collision integral [15]. 



In this paper, we will use the lattice Boltzmann method (LBM) [4, 37], which can be 
viewed as a special discretization of the Boltzmann equation [16]. In an LBM method, the 
distribution function f(x,v,t) is discretized on a space-time lattice with grid spacing dx in 
space and dt in time. Only a discrete number of velocities are considered, which correspond 
to a movement over an integer number of lattice points during one time-step, 

dx 
^i = Ci-^, Ci = -q,-q + l,...,q-l,q. 

For ease of exposition, we restrict ourselves to the case g = 1, which gives only three speeds 
(the so-called D1Q3 model). 

For reaction-diffusion systems, we can write the evolution law for fi{x, t) = f{x, Vi, t) as 

fi{x + Cidx,t + dt) = {l-uj)fi{x,t)-Lof-%x,t)+Ri{x,t), i = -1,0,1. (14) 

The right-hand side approximates the collision operator, and is composed of a reaction term 
Ri(x,t) and a BGK relaxation to the local equilibrium. After collision, the post-collision 
values are propagated to a neighboring lattice site, which corresponds to free flight. 

For the lattice Boltzmann discretization, we can define the (non-dimensional) moments 
of the distribution function as 

1 1 ^1 

P{x,t) = ^ fiix,t), 4>{x,t) = ^ Cifi{x,t), Cix,t) = 2 X] '^*^/«(^'*)- (1^) 

i=—l i=—l i=—l 

It can easily be verified [5, 29] that, under suitable smoothness assumptions, the system is 
well approximated by a macroscopic reaction-diffusion equation 

dtp{x,t) = d^{Dd^p{x,t)) +r{p{x,t)), (16) 

which can again be derived from the LBM equation using a Chapman-Enskog expansion 
[4, 42] using 

^= l^SDdt/dx^ ' ^.(^,i) = yr(p(x,i)), f^\x,t)=p{x,t)/3. (17) 
To view these models in the equation- free context, we define 

ti(x, t) = u{x, t) = {/i(x, t)} ■=_! , t/(X, t) = U{x, t) = p{x, t). 

3.2 Model problems and traveling waves 

In this paper, we propose a numerical method to compute coarse-grained traveling wave 
solutions of lattice Boltzmann models of the form (14). Traveling waves are solutions of the 
form 

U{x,t) = U{x-ct) = U{C), lirn U{0 = U^, (18) 



where C, = x — ct. They appear as steady states in a moving frame {C,-,t). 

In this work, we consider travehng fronts, for which U^ 7^ ?7^. For these sohitions, 
the evolution of the density is not sufficiently well described by the approximate reaction- 
diffusion PDE (16), due to a lack of smoothness of the density in the reaction front. Never- 
theless, we also observe in this case that the higher order moments (momentum and energy) 
become (more complicated) functionals of the density on fast time-scales. As a consequence, 
a PDE for the density should still exist. 

We consider two model problems: a model which is derived from the Fisher PDE, and a 
lattice Boltzmann model for ionization. Both model problems exhibit traveling fronts with 
arbitrary speed €> c* , where c* is called the critical speed. 

Model problem 3.1 (Fisher model). We consider the reaction-diffusion lattice Boltz- 
mann model, 

fi{x + Cidx, t + dt) = (1 - u)fi{x, t) - iof-\x, t) + Ri{x, t), (19) 

where i = —1,0, 1, with 

Ri{x,t) = ^p{x,t){l-p{x,t)). (20) 

Following the reasoning of the previous section, we define the approximate macroscopic 
PDE, 

dtp{x, t) = d^{Dd^p{x, t)) + p{x, t){l - p{x, t)), (21) 

which was originally proposed by Fisher as a model for the spread of advantageous genes 
[11]. This equation appears in a range of physical and biological applications exhibiting 
waves, see e.g. [28], and supports traveling fronts of the form p{x,t) = p{C), with 

lim p{0 = l, lim piO=0. 

C,~>~oo ^— >oo 

This model problem can be put in the equation-free framework by identifying 

^(x, t) = u{x, t) = {fi{x, t)}|=_i , t/(X, t) = U{x, t) = p{x, t). 

The Fisher equation is arguably the simplest model problem that exhibits traveling 
waves of arbitrary speed. It is known for the corresponding PDE that traveling waves exist 
with speed c > c* = 2\/D. Numerically, traveling waves exists for all c. However, the 
critical speed c* is the lowest speed for which the traveling waves are uniformly positive. 
The coarse-grained traveling wave with minimal speed c* = 2\/D is shown in figure 1. 

Model problem 3.2 (Planar streamer fronts). Streamers are sharp, non-linear waves 
of electrons that propagate through gases in the presence of strong electric fields. The strong 
field accelerates the electrons that cause ionization reactions during the collisions with the 
neutral gas particles. This ionization reaction creates additional electrons that are, again, 
accelerated. This results in an avalanche at the tip of the wave front. 




Figure 1: Coarse-grained traveling wave solution of the Fisher lattice Boltzmann model (19)-(20). 



We consider a lattice Boltzmann model for the distributions fi(x,t) of electrons, which 
is coupled to a PDE that governs the evolution of the electric field through the electron 
density p{x,t) = ^^ fi{x,t). The coupled equations are 



fi{x + Cidx,t + dt) 
dtE{x,t) 



(1 - Lo)fiix, t) - ujfl\x, t) - E{x, t) Zj V^JfJ{x, t) + R,{x, t), 
-p{x, t)E{x, t) - Dd^p{x, t) 



(22) 
in which the electric field E{x, t) appears as an external force. In the original Boltzmann 
equation (8), this force appears as E{x,t)dvf{x,v,t); in the lattice Boltzmann equation, 
this external force is discretized as E{x,t) ^ ■ Vijfj, as proposed by Luo [23]. 

In general, the reaction terms Ri{x,t) should be modeled on a microscopic level. How- 
ever, due to reasons of computational complexity, most analysis of streamers is based on 
the Townsend approximation of the microscopic ionization reactions that appear at the tip 
of the front [9]. The reaction terms are then given by Ri{x,t) = dtp{x,t)g{E{x,t))/3, with 
g{E) = \E\exp{-l/\E\). 

In this paper, we will also use this Townsend approximation to analyze the performance 
of our numerical method. In a forthcoming publication, the method presented here will be 
used to analyze the traveling waves of a five-speed lattice Boltzmann model, which is based 
on a more realistic set of microscopic interactions [43] . 

Again, we can derive an approximate PDE model, 

r dtp{x, t) = d^{p{x, t)E{x, t)) + Ddlp{x, t) + p{x, t)g{E{x, t)), 
\ dtE{x,t) = -p{x,t)E{x,t) - Dd^pix,t), 

This coupled equation exhibits traveling wave solutions with arbitrary speed c > c* 
limT. 



(23) 



\E{x)\ + 2^Dg{E{x)), see figure 2. 
This model problem can be put in the equation-free framework by considering 

n(x, t) = u{x, t) = {{f,{x, t)]]^_^,E{x, t)] , ?7(X, t) = U{x, t) = {p{x, t),E{x, t)} . 
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Figure 2: Coarse-grained traveling wave solution of the lattice Boltzmann model for planar streamer 
fronts (22). We show both the density p{x,t) (left) and the electrical field E{x,t) (right). 

4 The fixed point problem 

To define the fixed point problem for the traveling wave, we first construct a coarse-grained 
time-stepper for the lattice Boltzmann model of the form 

U{x,t + 6t) = S{U{x,t),6t), (24) 

where the coarse-grained time-step 6t = kdt, see equation (4). In section 4.1, we describe 
the details of the lifting operator. 

Traveling waves appear as a one-parameter family of solutions: together with U*{C), 
any translate U*{C + ip), for an arbitrary but fixed 99, is also a traveling wave. We therefore 
add a phase condition to remove this indeterminacy. The specific implementation of the 
shift-back operator and the phase condition are discussed in section 4.2. 

Note that we consider model problems for which a traveling wave exists for arbitrary 
c > c*. It is of physical interest to study how the critical speed c* depends on the system 
parameters. However, in this paper we confine ourselves to the computation of a traveling 
wave for a fixed speed c. A detailed analysis of the physical properties of the ionization 
model can be found in [43]. 

4.1 Lifting for lattice Boltzmann models 

As already outlined in section 2, a coarse-grained time-stepper is constructed as a lift- 
simulate-restrict procedure. While the restriction step is well defined by equation (15), the 
lifting step can be defined in multiple ways. We discuss three possibilities. 

Weighted lifting. A first possibility is to simply initialize the distributions as 

fi{x,t) =Wip{x,t), '^Wi = l, (25) 
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where Wi = 1/3 is an obvious choice, since this corresponds to the local diffusive equilibrium 
as defined in equation (17). As a consequence, the higher order moments are initialized as 

(A(x,t) = 0, ^{x,t) = pix,t)/3. 

This very rough approximation of the higher order moments introduces a lifting error, and 
one relies on the separation of time-scales between low order and high order moments to 
heal this error quickly [38] . In [40] , it is shown how this initialization can produce undesired 
artifacts for lattice Boltzmann models when 5t is small. 

Slaving relations. From the assumption that a macroscopic model for p{x, t) exists, it 
follows that the higher order moments (j){x,t) and ^{x,t) can be written as functionals of 
p{x,t). For lattice Boltzmann models of the form (14) with (17), these so-called slaving 
relations can be written down analytically as an asymptotic expansion in l/w. Up to third 
order, we have [41], 

, , , 2dx „ , , dxdt ( uj l\/o /, NN o2/ N\ 

nx^t) = -^—d^p{x,t) + — 2- ^+ Q {dxr{p{x,t))-d^tp{x,t)) , 

i{x, t) = -p{x, t)--—{r {p{x, t)) - dtp{x, t)) . 
6 buj 

These expansions can alternatively be written down in terms of p{x, t) and its spatial deriva- 
tives only, by making use of (16). The weighted lifting scheme coincides with the zeroth 
order term of the slaving relations. 

Constrained runs scheme. Although ideally one would like to use the slaving relations 
to initialize the lattice Boltzmann time-stepper, this is not always possible. The analytical 
derivation is cumbersome, and might even be impossible when the lattice Boltzmann model 
is coupled to a PDF, as in model problem 3.2. Moreover, the number of higher order 
moments increases when a more detailed discretization of the velocity space is used. Finally, 
the analytic derivation of the slaving relations is only tractable if only a small number of 
terms of the asymptotic expansion are needed. However, the macroscopic equation becomes 
invalid precisely when a large number of terms is needed. 

As an alternative, one can use the constrained runs scheme [12, 13] to approximate the 
full microscopic state corresponding to a set of predescribed macroscopic variables. The idea 
is to perform a number of lattice Boltzmann time-steps using equation (19), where after each 
time-step the density is reset to the initial density. During this iteration the microscopic 
variables converge towards their relaxed values; correspondingly, the higher order moments 
have then approximately reached the slaved state. As such, the constrained runs scheme is a 
fixed point iteration for the higher order moments of the microscopic model. We reproduce 
the schematic representation that was given in [41] (figure 3). 



11 



{<^(0),^(0)| __ 



{</.(p(o)),e(p(°))}t 




■ ■ . . approximate 
slow 



Figure 3: Sketch of the evolution of the constrained runs scheme. We plot the higher order moments 
(f> and ^ are plotted as a function of the density p. After each lattice Boltzmann step, p is reset to 
the given p^^\ The successive iterates converge to {(p*,^*}, which is an approximation to the slaved 
state {(/>(p^°^),C(p'"')} as given by equation (26). Figure reproduced from [41]. 

The constrained runs scheme can readily be applied to any lattice Boltzmann model, and 
defines a lifting operator that initializes the microscopic variables very close to their relaxed 
values. We refer to [41] for a detailed convergence and error analysis for reaction-diffusion 
systems. 



4.2 Shift-back procedure and phase condition 

Next, we discuss our specific implementation of the shift-back operator a^p. By noting that 
a shift-back over a distance y? is equivalent with a time evolution of the transport equation 

dtU{x,t)-d^u{x,t) = o, 

over a time 93, we obtain 

a^ : U{x, t) ^ a^{U{x, t)) = U{x, t) + ipd^U{x, t), (27) 

which is a valid approximation for shift-backs over short distances. 
At this point, we have constructed a non-linear system 

U{C)-ac5t{S{U{O,St))=0, (28) 

Unless we add a phase condition, this system is singular. When defining the Jacobian of 
the coarse-grained time-stepper as 



L{UiC),St) 



d[a,st{S{U{0,St))] 



duic) 

this singularity appears as a zero eigenvalue of I—L{U*{C),6t), with eigenfunction dU*{C)/dC- 
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To compensate for the extra phase condition, we add a regularization parameter a as an 
additional unknown, similar to what is done in [33] for Hamiltonian systems. The resulting 
non-linear system is 



G(^,a) 



V -acSt[S{\J,^t))^aAcP = 0, 
p(C/) = 0, 



(29) 



where \J denotes the space discretization of t/(C) and p(C/) is the phase condition defined 
below. This problem is well-posed for the unknowns (C/, a) with a locally unique solution 
(t7*,a* = 0). 

Here, we use the phase condition 

/"CiV-l 

vm = / i7(C)dcC/.e/(C)dC, 

-'Co 

which minimizes phase shift with respect to the reference solution Uref{C)i &s is done in 
publicly available software for bifurcation analysis, such as AUTO (for ordinary differential 
equations) [6, 7, 8] or DDE-BIFTOOL [10] (for delay differential equations). 

5 Preconditioned Newton-GMRES 



We solve the non-linear equation (29) for the traveling wave solution. First, we discretize 
C7(C) on the lattice Boltzmann grid [xo,xi = xq + dx, . . . ,xn-i], and provide the time- 
stepper with homogeneous Neumann boundary conditions. The spatial derivative in (27) is 
discretized using central differences. 

We solve the resulting non-linear system (29) with a Newton-Raphson procedure 



a(fc+i) 






(30) 



where dU^ ' and da^^ are the corrections calculated each iteration by solving linear systems 
of the form 



A{u'-''\6t) 






I-L{U^^\5t) d^tJ^^'^ 
dijp {tj^^'^) 



-G(C/(^^ 



a 



(kh 



(31) 



where A(C/'^,(^t) denotes the linearization of G{U,a) around the point {U^'^\a^^'). The 
right hand side is the residual at the same point. 

We do not have an explicit formula for L[U^^' ,6t), since it involves the Jacobian of 
the coarse-grained time-stepper with shift-back. However, we can estimate matrix-vector 
products as 

a,st {S{U + ev, 6t)) - a,st {SiU, 6t)) 



{l-L{U,6t)) 



V K, V 



(32) 
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Therefore, we solve the hnear system (31) using a Krylov method, such as GMRES [32]. 

The convergence rate of GMRES depends sensitively on the spectral properties of the 
system matrix in equation (31). For GMRES to converge rapidly, the eigenvalues should 
be clustered, e.g. around one [31]. It can be checked that the bordering row and column 



transform the zero eigenvalue of the Jacobian of (28) into itW dgp ■ d(^U, while leaving the 
other eigenvalues unaltered [25]. The eigenvalues //fe of 1 — L{U, 6t) are of the form 

/ifc = 1 -exp(Afc(5i), 

where A^ ~ 0{—k'^) for A: S> 1. For k small, Aj; ~ 0. When the system possesses a low- 
dimensional inertial manifold, one can choose 6t large, such that the spectrum becomes a 
compact perturbation of the unit matrix [18]. In this case, GMRES is known to converge 
very quickly; this, however, at the cost of a long simulation time for each matrix-vector 
product. If 6t is chosen to be small, the eigenvalues ^u^ range from to 1. In that case, 
preconditioning will be necessary. 

We define a preconditioning matrix M{U,6t), and replace the linear system (31) with 



M{U^^\6t)-^A{U^''\5t) 






-M{U^^\6t)-^G{U^''\a^''^] 



(33) 



Ideally, M(U, 5t) is both a good approximation of the system matrix, as well as a matrix for 
which an efficient (direct) solver is available. We will use a time-stepper for the approximate 
macroscopic equation (12) to construct M{U,6t). Consider the macroscopic equation (12) 
in the moving frame (Ci*)) 



dtUic, t) = F (u{c, t), d^u{c, t), . . . , dpic, t)) + cd^uic, t), 

= FiUiC,t)). 



(34) 



We construct a central finite difference/backward Euler time-stepper for (34) as 

U{t + dt) = S {U{t), St) = [I- 6tJ (C/(t))]~^ U{t), (35) 

where again U{t) is the space discretization of U{C, t) on the lattice Boltzmann grid, and 

dPiUit)) 



J{U{t)) 



dU 



We then define 



M{U,6t) 



I -{I- 6tJ{U)) dt^U 
dopitJ) 

with which we will solve linear systems of the form 



M{U,6t) 



dU 
da 



Ju 



(36) 



(37) 
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The first A^ equations can be simplified through some algebraic manipulation, 



I-{l-6tJ{U) dU + da-d(;U = bfj, 

dU ={I- dtJ{Uy^ dtJ + bu-da- d^f/, 

dU ={I- 5tJ{uy^ dU + bu, 
{I - 5tJ{U)) dU = dU+{l- 6tJ{U)) b(j, 
6tJ{U)dU = -{!- 6tJ{U)) bu, 
which leads to the equivalent linear system 

(I 



5tJ{U) 
djyp{tj) 



5tJ{U)) d^U 




dU 
da 



I - 5tJ{U)) bfj 
ba 



(38) 



in which the system matrix is the sum of a band and a rank one matrix. There exist efficient 
(order A^) direct solvers for linear systems of this type, see e.g. [3]. 

6 Numerical results 

6.1 Fisher equation 

We consider the model problem 3.1 on the domain [0, L] with L = 10 and with diffusion 
coefficient a = 0.1, and grid parameters dx = 2.5 ■ 10~^ and dt = 1 ■ 10"'^. For these 
parameters A^ = 400 and uj ~ 1.35. We intend to find the traveling wave with minimal 
speed c* = 2\AD, which is depicted in figure 1. The initial guess U^^' for the Newton 
process is obtained by a time integration of the Fisher PDE (21) over a time interval of size 
9000dt with initial condition 



U{x, 0) = p{x, 0) = exp(-5x''), X G [-L, L], 



(39) 



from which we take the part where x > 0. 

In order to choose an adequate value for the coarse-grained time-step 5t, we first il- 
lustrate the effects of the lifting operator. We then investigate the spectrum of the linear 
system that we need to solve, along with the effect of the preconditioner, and we conclude 
by showing the convergence of the Newton process. 

6.1.1 Initialization 

To understand the details of the initialization, we compare the exact evolution of the Fisher 
lattice Boltzmann model with the evolution after re-initialization. We perform an initial 
simulation of 100 steps with the lattice Boltzmann model, starting from the initial condition 
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Figure 4: Evolution of the initialization error for the Fisher problem for different lifting procedures. 
We show the difference between the distribution function of the re-initialized and the original lattice 
Boltzmann simulation. 

U^^' , which is lifted to distribution functions using the first order slaving relations. We 
extract the density at t = 80dt using equation (15). With this density, we initialize a 
second lattice Boltzmann simulation using the three procedures described in section 4.1. 
For the slaving relations, we only use the first order approximation. We then compare the 
evolution of the distribution functions with the evolution of the original simulation as we 
continue to evolve both. The result is somewhat unexpected. 

In figure 4, we plot the error between the original solution and the re-initialized simu- 
lation in the first 20 time-steps after re-initialization. We see that, in line with the results 
reported in [41], both the constrained runs scheme and initialization using first order slav- 
ing relations lead to a reduced lifting error, compared to weighted lifting. Unexpectedly, 
however, we observe that after initialization, all re-initialized simulations show an initial 
convergence towards the exact solution with a convergence rate of |1 — w|. This convergence 
stagnates after approximately 10 time-steps for this choice of to. This behaviour is not 
completely understood yet. As is to be expected, the smallest error is obtained by fully 
converged constrained runs. 

These results suggest that the simulation phase of the coarse-grained time-stepper, which 
involves the evolution of the lattice Boltzmann model from t* to t* + 5t requires a coarse- 
grained time-step 6t that is at least 10 times the microscopic time-step dt to eliminate 
initial transients caused by the lifting for this choice of uj. If a shorter time is choosen, 
undesired artifacts will show up in the spectrum of the coarse-grained time stepper. In our 
simulations, we choose 6t = 15dt = 1.5 • 10"^, and we initialize using the constrained runs 
scheme with 10 steps. 
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Figure 5: Fisher traveling wave problem. Left: Spectrum of the system matrix of equation (31). 
Right: spectrum of the preconditioning matrix (36). 

6.1.2 Performance of the preconditioner 

We now compute the spectrum of the matrix A{lJ^^\^t\ using the parameters defined 
above. To this end, we construct the matrix Aijj^^' ^bi) explicitly by computing the matrix- 
vector products with all coordinate vectors ej = [0, . . . , 0, 1, 0, . . . , 0]-^. The matrix-vector 
products are estimated using (32) with e = 1 • 10^^. The results are shown in figure 5 (left). 
We compare with the spectrum of the matrix M(\J^^\ 6t) shown in figure 5(right). We see 
good agreement between the two spectra. The differences are in the imaginary part of the 
eigenvalues, and in the very fast modes, which show up around 1 since we are computing the 
spectrum of the fixed point iteration. Also remark how the zero eigenvalue, corresponding 
to the singularity of the fixed point equation (28), is split into two isolated eigenvalues by 
the addition of the phase condition and the artificial unknown a. 

The spectrum of the matrix M{U^^' ,6t)~^A{U^^' ,6t) is shown in figure 6(left). We see 
that this spectrum is nicely clustered around 1. As a result, we obtain very fast GMRES 
convergence. Figure 6 shows the error as a function of the number of GMRES iterations. 
We clearly see the fast linear convergence, as predicted by the theory [31]. 



6.1.3 Convergence of the Newton process 

We are now ready to show the convergence of the Newton process. To this end, we take as 
an initial guess U^^' a smoothed step function of the form 

1 



U(0) 



(40) 



exp{2{x-L/2)) + V 

with L = 10, the length of the interval. The convergence of the Newton process is shown 
in figure 7. We notice that the convergence is linear, and the convergence factor changes 
after approximately 4 iterations. In the process, we have converged to the traveling wave 
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Figure 6: Fisher traveling wave problem. Left: spectrum of the product of the inverse of the 
preconditioner and the system matrix. Right: error of the hncar system as a function of the number 
of GMRES iterations. 

solution shown in figure 1, and a* = —2.24 • lO"'^. The fact that a* is not equal to zero can 
be shown to be an artifact of the space discretization, combined with the truncation to a 
finite domain. Indeed, a* becomes smaller when we refine the lattice Boltzmann grid (and 
decrease the time-step to keep uj fixed). We can explain this as follows: the conclusion that 
a* should equal zero follows from the assumption that the extra column that is added in the 
Jacobian, d^C/, is the eigenvector associated with the zero eigenvalue of the fixed point map 
(28). However, due to the space discretization, d(^U will only be an approximation to this 
eigenvector. Therefore, we are using a quasi-Newton method (hence the linear convergence), 
and the solution value a* will be non-zero but small. The convergence factor changes when 
the traveling wave solution has been computed to machine precision. At that point, only 
the error in a shows up in the residual, indicating that the non-zero value of a* is indeed 
the main cause for the non-quadratic convergence. 



6.2 Planar streamer front 

As a second example, we numerically study model problem 3.2 for planar ionization fronts 
on the domain [0, 130] with grid parameter dx = 1 • 10^^ and dt = 1 ■ 10"'^. We choose 
the diffusion coefficient to be a = 1 and the electric field E^ = —1 (i.e. a constant, large 
electric field on the far right). For these parameters, the corresponding lattice Boltzmann 
relaxation parameter is u; ~ 1.538, and N = 1300, which implies that the unknown solution 
has 2601 unknowns. We intend to compute the traveling wave solution with the minimal 
speed c* = \E^\ + 2\/ Dg{E^), which is depicted in figure 2. The initial guess U^^' for the 
Newton process is obtained by a time integration of the ionization PDE (23) over a time 
interval of size 4000dt with initial condition 



U{x,0) = pix,0) = exp{-15x^),x€ [0,130]. 



(41) 
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Figure 7: Fisher traveling wave problem. Convergence of the Newton process as a function of 
iteration number. 
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Figure 8: Evolution of the initialization error for the planar streamer model with different lifting 
procedures. We show the difference between the distribution function of the re- initialized and the 
original lattice Boltzmann simulation. 

6.2.1 Initialization 

Again, we study the properties of the initialization by comparing the exact evolution of 
a lattice Boltzmann simulation with the evolution after re-initialization. We performed a 
simulation of 100 steps, starting from U^^' , with the lattice Boltzmann model (22), which 
is initialized using the first order slaving relations. Again, we extract the density at time 
t = 80dt, and initialize a second lattice Boltzmann simulation using the three procedures 
described in section 4.1. We observe the same behaviour as for the Fisher model. In 
figure 8, we show the evolution of the initialization error during the first 20 steps after 
re-initialization. We again see an initial convergence towards the exact solution, which 
stagnates after 10 to 20 time-steps. Depending on the accuracy of the lifting step, we get a 
better correspondence with the correct distribution functions. 
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Figure 9: The planar streamer front problem. Left: Spectrum of the system matrix of equation 
(31). Right: spectrum of the preconditioncr (36). 

Based on these observations, we choose bt = 15df = 1.5 • 10^^ in our simulations, and 
we initialize using the constrained runs scheme with 10 steps. 

6.2.2 Performance of the preconditioner 

In figure 9, we show the spectrum of the matrix A{tj^^' ,6t), using the parameters defined 
above, and the spectrum of the corresponding preconditioning matrix M{tj^^' ,5t). Again, 
we see good qualitative agreement between the two spectra, but the difi'erences are more 
pronounced than for the Fisher equation. 

The spectrum of the matrix M{U^^\5t)~^ A{tj^^\5t) is shown in figure lO(left). Al- 
though most eigenvalues are in a small cluster around 1, a few eigenvalues are scattered on 
the real axis. These additional eigenvalues result in a temporary stagnation of the GMRES 
convergence, as shown in figure 10 (right). Between these stagnations, linear convergence 
is observed. We note that, even with the temporary stagnation, the solution is computed 
in 30 to 40 iterations, while the number of unknowns is 2601. 

6.2.3 Convergence of the Newton process 

We proceed to show the convergence of the Newton process, starting from the initial guess 
U^^' . The results are shown in figure 11. We plot the norm of the residual and the correction 
{dU^ ' ,a^ ') after each iteration. For this problem, we get quadratic convergence towards 
the solutions shown in figure 2, and a* = 0. Thus, in this case, the effect of the discretization 
is much less pronounced, or even absent. 
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Figure 10: Planar streamer front problem. Left: spectrum of the preconditioned system matrix 
M~^A (equation (33)) for the first Newton step. Right: the convergence of the preconditioned 
GMRES procedure as a function of the number of iterations. 
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Figure 11: Planar streamer front problem. Convergence of the Newton process as a function of 
iteration number. 
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7 Discussion and conclusions 

In this article, we showed how one can use a coarse-grained time-stepper to compute trav- 
ehng wave solutions of lattice Boltzmann models. In a co-moving frame, emulated by 
performing a shift-back operation after each coarse-grained time-step, the traveling wave ap- 
pears as a steady state, which is computed using a Jacobian-free Newton-GMRES method. 
The method uses repeated calls to the coarse-grained time-stepper to estimate the required 
matrix- vector products. For efficiency reasons, we limited the size of the coarse-grained 
time-step 6t. The real part of the spectrum of the resulting Jacobian ranges from to 1, 
which results in slow convergence of Krylov methods. Therefore, we accelerated convergence 
of the GMRES procedure by introducing a preconditioner that is based on an approximate 
macroscopic model, which is derived from the lattice Boltzmann model using a Chapman- 
Enskog expansion. We illustrated the approach on a lattice Boltzmann model for the Fisher 
equation and on a model for an ionization wave. 

The total cost of finding a solution is determined by the required number of microscopic 
time-steps, which depends on the number of GMRES iterations and on the number of 
microscopic time-steps per matrix-vector product. For each matrix-vector product, we use 
the microscopic time-stepper both during the simulation and the lifting step when the 
constrained runs scheme is used. The required number of calls to the microscopic time- 
stepper is determined by the microscopic relaxation parameter lo, since the convergence 
rate of the distribution functions towards their correctly slaved values is given by 1 1 — u;| . In 
our examples, approximately 10 lattice Boltzmann time-steps were needed for each lifting 
operation and an additional 10 to 15 for each microscopic simulation. 

The convergence rate of the GMRES method is closely related to the quality of the 
preconditioner. The preconditioner performs well when the approximate macroscopic model 
is a good approximation of the lattice Boltzmann model. Our experiments indicate that 
the quality of the preconditioner is not related to the number of lattice points N, which 
implies that the number of GMRES iterations is approximately constant with varying A^. 

Since the number of time-steps only depends on the relaxation parameter uj and the 
spectrum of the preconditioned time-stepper, we conclude that the convergence rate is 
independent of the number of variables in the problem. As a consequence, the algorithm 
scales with the cost of taking a single lattice Boltzmann time-step and the cost of a direct 
solve with the preconditioning matrix. 

In our numerical examples we found that, on average, 2000 time-steps are required 
to converge to the traveling wave solution. This is still expensive if we compare with 
straightforward time integration. E.g. for the Fisher model, a localized initial state evolves 
into a steady traveling wave within 9000 time-steps. However, direct time integration looses 
its appeal in situations where the microscopic time-step is so small that it is computationally 
infeasible to reach the time horizon where the traveling wave becomes steady. Furthermore, 
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the total number of time-steps will be much smaller in the context of continuation, where 
the behaviour of the solutions is studied as a function of one or more varying parameters. 
In this setting, the solution for nearby parameter values provides an accurate initial guess, 
hereby lowering the required number of Newton steps. Also, the Newton-GMRES method 
allows us to find unstable traveling waves. These will never be found by simulation, since 
any perturbation will grow, and will destroy the traveling wave. 

Note that preconditioning with a roughly approximate macroscopic PDE model can 
also be helpful to study the linear stability properties of the traveling wave. Indeed, the 
relevant rightmost eigenvalues of the Jacobian matrix in the solution can be computed by 
an Arnoldi iteration that requires the solution of a linear system in each iteration [21]. This 
linear system can again be preconditioned with the techniques proposed in this article. 

In principle, the Newton-GMRES method could be applied to the lattice Boltzmann 
model directly. However, it is much harder to find an appropriate preconditioner in that 
setting. Indeed, one cannot expect a preconditioner based on the approximate PDE to 
work well, since the fast time-scales, which describe the slaving of the higher order mo- 
ments, will not be correctly accounted for. We also emphasize that, from the coarse-grained 
solution, the detailed solution can easily be reconstructed using the constrained runs lifting 
scheme. Our method generalizes readily to lattice Boltzmann models with a more detailed 
description of the velocity space. A particular five-speed model, in which the Townsend 
approximation for the reaction term is replaced by a more realistic set of microscopic inter- 
action rules, is currently being investigated [43]. 
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